20180716_GB_SH1_d13Call_presentation

Load data

#> Info: preparing to read 1 data file(s)...
#> Info: reading file 1/1 'Garrett_Adducts.cf.rda' with '.rda' reader...
#> Info: loaded data for 71 data files from R Data Archive - checking loaded files for content consistency...
#> Info: encountered 6 problems in total.
#> # A tibble: 6 x 4
#>   file_id                  type  func           details                   
#>   <chr>                    <chr> <chr>          <chr>                     
#> 1 587__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 2 587__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
#> 3 586__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 4 586__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…
#> 5 585__Linearity CO2__.dxf error extract_dxf_r… cannot identify measured …
#> 6 585__Linearity CO2__.dxf error extract_isoda… `peak`, `start`, `rt`, `e…

Vendor data table

#> Info: aggregating vendor data table without units from 71 data file(s), including file info 'c(ox = `Seed Oxidation`, everything())'

Map peaks

Chromatograms

Retrieve metadata

Map peaks

# add metadata
data_w_metadata <- vdt %>% 
  iso_add_metadata(metadata_samples, match_by = c(`Identifier 1`)) 
#> Info: metadata added to 1827 data rows, 40 left without metadata:
#> - 25 metadata entries were mapped to 1827 data entries via column 'Identifier 1'

# missing metadata
data_w_metadata %>% 
  iso_get_missing_metadata(select = c(Analysis, `Identifier 1`, map_id))
#> Info: fetching data entries that are missing metadata

# map peaks
data_w_peaks_all <- data_w_metadata %>% 
  # only continue with records that have metadata and are flagged for processing
  iso_remove_missing_metadata() %>% 
  filter(process == "yes") %>% 
  # map peaks
  iso_map_peaks(metadata_peak_maps, file_id = file_id, rt = Rt, rt_start = Start, rt_end = End)
#> Info: removing 40 of 1867 data entries because of missing metadata
#> 

# problematic peaks
prob <- data_w_peaks_all %>% 
  iso_get_problematic_peaks(select = c(file_id, Analysis, compound, Start, Rt, End))
#> Info: fetching 596 of 2044 data table entries with problematic peak identifications (unidentified, missing or ambiguous)

# focus on non problematic peaks only?
data_w_peaks <- data_w_peaks_all %>% iso_remove_problematic_peaks()
#> Info: removing 596 of 2044 data table entries because of problematic peak identifications (unidentified, missing or ambiguous)
data_w_peaks

Reference peaks

data_w_peaks %>% 
  iso_plot_ref_peaks(
    is_ref_condition = is_ref_peak == "yes", 
    ratio = c(`R 13C/12C`, `R 18O/16O`),
    group_id = `Analysis`,
    within_group = TRUE
  ) 

Data processing

Focus on analytes

Focus on the analytes and calculate a few summary parameters we want to use later.

Evaluate Standards

Known isotope values

Add known isotope values for standards.

#> Info: matching standards by 'type' and 'compound' - added 15 standard entries to 458 out of 762 rows

Single analysis calibration (for QA)

Determine calibrations fits for individual standard analyses.

#> Info: preparing data for calibration by grouping based on 'Analysis'
#> Info: generating 'd13C' calibration based on 1 model ('lm(`d 13C/12C` ~ true_d13C)') for 31 data group(s) with standards filter 'is_standard'. Storing residuals in new column 'd13C_resid'.
#> Info: there were no problematic calibrations
#> Warning in `[<-.data.frame`(`*tmp*`, is_list, value = structure(list(`3` =
#> "<>", : replacement element 1 has 1 row to replace 0 rows
#> Warning in `[<-.data.frame`(`*tmp*`, is_list, value = structure(list(`3` =
#> "<>", : replacement element 2 has 1 row to replace 0 rows

Calibration coefficients

Look at calibration parameters (coefficients and summary) as well ad data ranges in overview table.

Visualize the calibration parameters

# unnest some of the data for plotting
data_for_plots <- 
  stds_w_calibs %>%
  iso_unnest_data(
    c(Oxidation = `Seed Oxidation`, 
      datetime = file_datetime,
      `Mean Amplitude` = ampl_sample_mean.mV)
  ) 
data_for_plots

# parameters vs time
data_for_plots %>%
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = datetime,
    color = `Injection Volume`,
    shape = Oxidation
  )

# parameter vs analysis
data_for_plots %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = paste(Analysis),
    color = `Injection Volume`,
    shape = Oxidation
  )

# parameters vs amplitude
data_for_plots %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C", 
    x = `Mean Amplitude`,
    color = `Injection Volume`
  ) 

# turn the last plot into an interactive one
ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Inspect Residuals

Look at residuals to see goodness of fit for the single analysis calibrations.

stds_w_calibs %>% 
  # pull out all data
  iso_unnest_data(select = everything()) %>% 
  filter(is_standard) %>% 
  # calculate deviation from means
  group_by(Analysis) %>% 
  mutate(
    `Var: residual d13C [permil]` = d13C_resid,
    `Var: area diff from mean [%]` = (`Intensity 44`/mean(`Intensity 44`) - 1) * 100,
    `Var: amplitude diff from mean [%]` = (`Ampl 44`/mean(`Ampl 44`) - 1) * 100
  ) %>%
  ungroup() %>% 
  # visualize
  iso_plot_data(
    x = compound,
    y = starts_with("Var"),
    group = Analysis,
    color = `Injection Volume`
  ) +
  # plot modifications
  labs(x = "")

ggplotly()

Gobal calibration

Test different calibration models

Run global calibrations across all standards.

#> Info: preparing data for calibration by nesting the entire dataset
#> Info: generating 'd13C' calibration based on 4 models ('delta_only = lm(`d 13C/12C` ~ true_d13C)', 'delta_and_ampl = lm(`d 13C/12C` ~ true_d13C + `Ampl 44`)', 'delta_cross_ampl = lm(`d 13C/12C` ~ true_d13C * `Ampl 44`)', 'delta_and_time = lm(`d 13C/12C` ~ true_d13C + file_datetime)') for 1 data group(s) with standards filter 'type == "standard" & is_standard'. Storing residuals in new column 'd13C_resid'.

Visualize the calibration parameters for the different models.

# coefficient summary table
data_w_calibs %>% iso_unnest_calibration_parameters("d13C")

# plot to look at coefficients and summary
data_w_calibs %>% 
  iso_plot_calibration_parameters(
    calibration = "d13C",
    x = d13C_calib,
    color = d13C_calib,
    # include the significance level to gauge which model parameters matter
    shape = signif
  )

Garrett thinks delta and aplitude and delta cross amplitude look better

data_w_calibs %>% 
  iso_unnest_data(select = everything()) %>% 
  rename(`Injection Volume` = `Identifier 2`) %>% 
  filter(type == "standard" & is_standard) %>% 
  # plot each standard analysis for each of the models to see the differences
  iso_plot_data(
    x = compound,
    y = d13C_resid,
    group = paste(Analysis, d13C_calib),
    color = d13C_calib
  ) 

Apply calibration

Apply the delta_and_ampl calibration.

#> Info: applying 'd13C' calibration to infer 'true_d13C' for 1 data group(s); storing resulting value in new column 'true_d13C_pred' and estimated error in new column 'true_d13C_pred_se'. This may take a moment... finished.

Inspect results

The comparison of the predicted value in the standards along with the standard error in the prediction (the inherent error of the global calibration itself, which should be treated as the analytical error) can give a better sense of how precise the data is.

Visualization

Keep in mind that the calibration provides not just the predicted value of a peak but also the standard error based on the calibration and an indicator whether a predicted data point was in the calibration range or not.

data_calibrated %>% 
  iso_plot_data(
    x = Analysis, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se,
    color = true_d13C_pred_in_range,
    shape = type,
    points = TRUE,
    lines = FALSE
  ) + 
  # additional features beyond the default plot
  facet_wrap(~compound, scales = "free_y") + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(color = "in calib. range")

Combine with other data for final plots

Remove bad chromtogram data

final_data <- data_calibrated %>%
  #filter out bad chromatograms - all are the lower amount injection, which were followed by a larger volume - keeping the larger volume injection data
  filter(Analysis != "597") %>% # 122 depth, OG005
  filter(Analysis != "593") %>% # 121.105, OG047
  filter(Analysis != "612") %>% # 121.89, OG040
  filter(Analysis != "613") %>% # 120.205, OG042-2
  filter(Analysis != "609") %>% # 119.55, OG158
  filter(Analysis != "595") %>% # 122.9, OG043
  filter(type == "sample") %>%
  filter(true_d13C_pred_in_range == TRUE) %>%
  select(
      # sample info
      depth = `Depth`,  `compound`, Identifier.1 = `Identifier 1`, `ampl_sample_mean.mV`,  type, Analysis, 
      # peak info
       Ampl.44 = `Ampl 44`,  true_d13C_pred , true_d13C_pred_se, true_d13C_pred_in_range,  d13C_resid
    ) %>%
  mutate(
    prep = "adduct"
  )

vis <- final_data %>% 
  filter(compound %in% c("C29", "C27", "C31", "C33"))%>%
  iso_plot_data(
    x = depth, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se,
    size = Ampl.44
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse() +
  coord_flip()
ggplotly(vis)

Low res alkane data

lowres <- read.csv("lowres_d13c.csv") %>% select( -X, -`measurement_info`, -`file_datetime`, -file_id) %>% mutate(prep = "non-adduct") 


sapply(final_data, class)
#>                   depth                compound            Identifier.1 
#>               "numeric"             "character"             "character" 
#>     ampl_sample_mean.mV                    type                Analysis 
#>               "numeric"             "character"             "character" 
#>                 Ampl.44          true_d13C_pred       true_d13C_pred_se 
#>               "numeric"               "numeric"               "numeric" 
#> true_d13C_pred_in_range              d13C_resid                    prep 
#>               "logical"               "numeric"             "character"
sapply(lowres, class)
#>                compound            Identifier.1                   depth 
#>                "factor"                "factor"               "numeric" 
#>     ampl_sample_mean.mV                    type                 Ampl.44 
#>               "numeric"                "factor"               "numeric" 
#>               true_d13C          true_d13C_pred       true_d13C_pred_se 
#>               "numeric"               "numeric"               "numeric" 
#> true_d13C_pred_in_range              d13C_resid                    prep 
#>               "logical"               "numeric"             "character"

final_data <- bind_rows(lowres, final_data)
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding factor and character vector,
#> coercing into character vector
#> Warning in bind_rows_(x, .id): binding character and factor vector,
#> coercing into character vector

Manipulate Alkane Data

Correct for biosynth fractionation

final_data <- final_data %>%
  filter(type == "sample") %>%
  select(-true_d13C, -d13C_resid) %>%
  filter(compound %in% c("phytane", "C17", "C19", "C21", "C27", "C29", "C31", "C33", "C35")) %>%
# net fractionation factor for inorganic -> lipid d13C (eps_net) calculated as = (inorganic -> biomass fractionation factor)  + (biomass -> lipid fractionation factor)

# minimum and max for marine algal biomass (first number in equation) from spread in data of algal biomass in Ohkouchi et al., 2015; OPTIONS FOR SHORTCHAIN: -4 per mil for phytane (Schouten et al., 1998; Sinninghe Damste et al., 2008), -8.4 per mil for cyano nC17 , -6.4 per mil for phytane (above takes into account frctionation in sulfur-bound and release per Sinninghe Damste?)(Sakata et al., 1997)
  
# lipid-leaf addition (second number in equation) by chain length from gymnosperm data in Diefendorf et al., 2011 [NOTE: average for gymnosperm is -2.5 per mil, that may be more appropriate]; minimum and maximum offset of leaf to atmosphere from C3 plant histogram in Ehleringer and Cerling, 2001 (in Taiz et al., 6th ed. Plant Physiology and Development)
  
  mutate(eps_net_min = case_when(.$compound %in% c("C17", "C19", "C21")  ~ (12+4) ,
                                 .$compound  == "phytane" ~ (12+4.4),
                                 .$compound  == "C27" ~ (22+2),  
                                 .$compound  == "C29" ~ (22+3),
                                 .$compound  == "C31" ~ (22+4),
                                 .$compound ==  "C33" ~ (22+3.5),
                                 .$compound ==  "C35" ~ (22+3)), 

         eps_net_max = case_when(.$compound %in% c("C17", "C19", "C21")  ~ (25+8) , 
                                 .$compound  == "phytane" ~ (25+6.4),
                                 .$compound  == "C27" ~ (32+2),  
                                 .$compound  == "C29" ~ (32+3),
                                 .$compound  == "C31" ~ (32+4), 
                                 .$compound ==  "C33" ~ (32+3.5),
                                 .$compound ==  "C35" ~ (32+3)),  
         
         d13C_pool_min = true_d13C_pred + eps_net_min,
         
         d13C_pool_max = true_d13C_pred + eps_net_max,
         
         d13C_pool_avg = true_d13C_pred + ((eps_net_max + eps_net_min) / 2) )

Calculate / estimate eps_p

final_data <- final_data %>%
    mutate(
         d13C_p = case_when(.$compound %in% c("C17", "C19", "C21")  ~ true_d13C_pred + 8.4,
                           .$compound == "phytane" ~ true_d13C_pred + 6.4 ) ,  # from Julio's unpublished text, original from Sakata et al., 1997 (6.4 in Sakata...)
         
         eps_p_lipid = 1000 * ((d13C_pool_avg + 1000)/ (d13C_p + 1000) -1), # from Freeman and Hayes 1992, but using avg estimated d13C of DIC from lipid rather than CaCO3 measurement. This allows moving pool value - important when pool itself is experiencing characterizable d13C changes thru time. 
         
         eps_f = case_when(.$compound %in% c("C17", "C19", "C21")  ~ 24,
                           .$compound == "phytane" ~ 25 ) ,    # from Julio's unpublished text, 22 for alkanes (cyano-plank intermediate) and 25 for phytane; Pagani 2002 says 25 for phytoplankton in general;  "Goericke et al. [ 1994] calculate the likely range of eps_f for marine phytoplankton to be 25-28" (from Brigidare et al., 1997). After Modelling (below), 24 for short chains gets us closest to phytane value (assuming 25 for phytane) - has implications for makeup of community, i.e., that cyanos are lower abundance than phytoplankton
         
         b = 170, # per Julio's unpublished text, can sequence from [118 to 262], or use constant 170. I think at least 170, potentially higher, given modelled need for high [PO4] advection to sustain OAE2
         
         CO2_aq_conc = b / (eps_f - eps_p_lipid) , # from Brigidare et al. 1997
         
         Ko = (2.418 * 10^-2), # solubility constant, from Weiss 1974 (tables II and III); dep. on T and Salinity, range from ~2 to ~7 over all regimes; @ ~35 C, range is from 2.1 to 2.3 across all salinities. Choosing for 34ppt (Hay et al., 1996) and 32C (mid range T from O'Brient et al., 2017)
         
         pCO2 = CO2_aq_conc / Ko  # from Julio's unpublished text , and Sarmiento and Gruber, Ocean Biogeochemical Dynamics, section 8.3 (p.327)
         )

Check calcultions immediately

final_data %>% filter(compound == "C21") %>% select( compound, eps_p_lipid, d13C_pool_avg, true_d13C_pred, CO2_aq_conc, pCO2) %>% kable()

compound	eps_p_lipid	d13C_pool_avg	true_d13C_pred	CO2_aq_conc	pCO2
C21	16.45642	-5.558408	-30.05841	22.53572	931.9983
C21	16.44228	-4.716870	-29.21687	22.49355	930.2542
C21	16.44119	-4.652097	-29.15210	22.49031	930.1204
C21	16.41665	-3.188455	-27.68846	22.41754	927.1108
C21	16.42721	-3.818942	-28.31894	22.44880	928.4037
C21	16.42893	-3.921605	-28.42161	22.45390	928.6147
C21	16.39894	-2.129421	-26.62942	22.36531	924.9508
C21	16.47220	-6.495370	-30.99537	22.58295	933.9514
C21	16.45913	-5.719790	-30.21979	22.54383	932.3338
C21	16.45878	-5.698920	-30.19892	22.54278	932.2904
C21	16.43330	-4.182080	-28.68208	22.46687	929.1508
C21	16.45231	-5.313985	-29.81398	22.52345	931.4907
C21	16.39212	-1.721012	-26.22101	22.34527	924.1218
C21	16.42510	-3.692964	-28.19296	22.44254	928.1449
C21	16.39335	-1.794682	-26.29468	22.34888	924.2712
C21	16.40776	-2.656710	-27.15671	22.39127	926.0244
C21	16.43223	-4.118406	-28.61841	22.46370	929.0197
C21	16.48668	-7.354039	-31.85404	22.62648	935.7518
C21	16.44045	-4.608185	-29.10818	22.48812	930.0296
C21	16.45807	-5.656454	-30.15645	22.54065	932.2021
C21	16.44026	-4.596741	-29.09674	22.48755	930.0060
C21	16.44202	-4.701641	-29.20164	22.49279	930.2227
C21	16.44346	-4.787331	-29.28733	22.49707	930.3999
C21	16.43946	-4.548936	-29.04894	22.48516	929.9073
C21	16.38667	-1.394377	-25.89438	22.32927	923.4603
C21	16.43286	-4.155748	-28.65575	22.46555	929.0966
C21	16.45481	-5.462670	-29.96267	22.53091	931.7994
C21	16.45315	-5.363880	-29.86388	22.52595	931.5943
C21	16.44810	-5.063625	-29.56362	22.51090	930.9717
C21	16.44502	-4.880004	-29.38000	22.50170	930.5916
C21	16.44379	-4.806987	-29.30699	22.49805	930.4405
C21	16.43948	-4.550064	-29.05006	22.48521	929.9096
C21	16.45443	-5.439936	-29.93994	22.52977	931.7522
C21	16.49080	-7.598096	-32.09810	22.63890	936.2654
C21	16.45631	-5.551668	-30.05167	22.53538	931.9843

Model role of measured and fixed variables in pCO2 output

data_try <- expand.grid(
b = seq(170, 262, by = 10), 
Ko = seq(2.1, 2.9, by = 0.2) * 10^-2,
eps_p = seq(16.5, 17.5, by = 1)
) %>% as_data_frame() %>%
  mutate(
    eps_f = case_when(eps_p ==16.5 ~ 24,
                      eps_p == 17.5 ~ 25),
    compound = case_when (eps_p == 16.5 ~ "C17,19,21",
                          eps_p == 17.5 ~ "phytane"),
    CO2_aq = b / (eps_f - eps_p),
    pCO2 = CO2_aq / Ko
)

# base plot
base_plot <- 
  ggplot() +
  aes(x = b, y = pCO2, linetype = factor(Ko)) + #removed color for eps_p because both eps_p lie on the same trend (becuase of eps_f assignment)
  geom_line() 

  
# plotting data 1
base_plot %+% data_try

Bring in complimentary data (Jones et al., 2018)

Inorganic geochem data

cisotope <- read_excel(file.path("metadata", "Appendix_Table1_geochemistry.xlsx")) %>% 
  #rename columns
  select(depth = `Abs. depth (m)` , d13c_org = `Average δ13Corg (‰ VPDB)` )

final_data <- full_join(cisotope, final_data, by = "depth")

Fit to floating timescale

timescale <- read_excel(file.path("metadata", "Appendix_Table2_SH1_agemodel.xlsx")) %>% 
  rename(depth = `Depth`) %>%
  arrange(depth)

# extrapolation
extrapolated_timescale <-
  data_frame(
    depth = final_data$depth %>% unique()
  ) %>% 
  filter(!is.na(depth)) %>% 
  mutate(
    in_range = depth <= max(timescale$depth) & depth >= min(timescale$depth),
    Date.Ma = map_dbl(depth, function(d) {
      fit_data <- mutate(timescale, closest = abs(d - depth)) %>% arrange(closest)
      m <- lm(Date.Ma ~ depth, fit_data[1:10,])
      p <- predict(m, data.frame(depth = d), se.fit = TRUE)
      return(p$fit)
    }) ,
    in_range = depth <= max(timescale$depth) & depth >= min(timescale$depth),
    ka = map_dbl(depth, function(d) {
      fit_data <- mutate(timescale, closest = abs(d - depth)) %>% arrange(closest)
      m <- lm(ka ~ depth, fit_data[1:10,])
      p <- predict(m, data.frame(depth = d), se.fit = TRUE)
      return(p$fit)
    }))  
    
ggplot() + 
  aes(depth, Date.Ma) + 
  geom_line(data = timescale) +
  geom_point(data = extrapolated_timescale, mapping = aes(color = in_range))

alldata_w_time <- left_join(final_data, extrapolated_timescale, by = "depth")

Estimate error for deposition

alldata_w_time <-  alldata_w_time %>%
  mutate(age_range = case_when(.$compound %in% c("C17", "C19", "C21", "phytane")  ~ (-0.00001) , #in MA (100 yrs)
                                 .$compound %in% c("C27", "C29", "C31", "C33", "C35") ~ (-.01)), #in MA (10000 yrs)
         age_dev = Date.Ma + age_range 
   )

Final plots

High res

alkane <- alldata_w_time %>% 
  filter(compound %in% c( "C29", "C31", "C33"), true_d13C_pred_in_range == TRUE, prep == "adduct")%>%
  iso_plot_data(
    x = Date.Ma, 
    y = true_d13C_pred, 
    #y_error = true_d13C_pred_se ,
    shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(94.6, 94.3))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(alkane)

bulkorg <-  subset(alldata_w_time, ka > -100 & ka < 120) %>%
  iso_plot_data(
    x = ka, 
    y = d13c_org
  ) + 
  # additional features beyond the default plot
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  #scale_y_reverse() +
  scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()


grid.arrange(bulkorg, alkane, ncol = 2)

Low res (all)

alkaneall <- alldata_w_time %>% 
  filter(compound %in% c("pristane", "phytane", "C19", "C21", "C31", "C33"), true_d13C_pred_in_range == TRUE)%>%
  iso_plot_data(
    x = Date.Ma, 
    y = true_d13C_pred, 
    y_error = true_d13C_pred_se ,
    color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse() +
  #scale_x_continuous()+
  scale_y_continuous() +
  coord_flip() 

ggplotly(alkaneall)

bulkorgall <-  subset(alldata_w_time) %>%
  iso_plot_data(
    x = ka, 
    y = d13c_org
  ) + 
  # additional features beyond the default plot
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  #scale_y_reverse() +
  #scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()


grid.arrange(bulkorgall, alkaneall, ncol = 2)

Pool composition estimates

pools <- alldata_w_time %>% 
  filter(compound %in% c("C19", "C21", "C31", "C33"))%>%
  ggplot()+
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_point(aes(x = Date.Ma , y = d13C_pool_avg)) +
  #geom_errorbar(aes(x = ka , ymin = d13C_pool_min, ymax = d13C_pool_max)) +
  #geom_errorbarh(aes(xmin = age_dev, xmax = ka , y = d13C_pool_avg)) +
  #geom_line(aes(x = age_dev, y = d13C_pool_avg, color = "red")) +
  #geom_line(aes(x = ka, y = d13C_pool_max)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "Inorganic pool d13C") +
  scale_x_reverse() +
  #scale_x_continuous()+
  scale_y_continuous() +
  #scale_x_continuous(limits = c(-100, 120)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()
ggplotly(pools)

Pool estimates with temporal offset estimates

pools_alttime <- alldata_w_time %>% 
  filter(Identifier.1 != "OG005_F1 alkanes_122m") %>% #removing duplicate for 94.45 Ma (non-adducted)
  filter(compound %in% c( "C19",  "C33"))%>%
  ggplot()+
  aes(x = Date.Ma , y = d13C_pool_avg, color = "red") +
  facet_grid(~compound, scales = "free", space = "free") + 
  geom_line() +
  #geom_errorbar(aes(x = ka , ymin = d13C_pool_min, ymax = d13C_pool_max)) +
  #geom_errorbarh(aes(xmin = age_dev, xmax = ka , y = d13C_pool_avg)) +
  #geom_line(aes(x = ka, y = d13C_pool_min)) +
  geom_line(aes(x = age_dev, y = d13C_pool_avg, color = "blue")) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "Inorganic pool d13C") +
  scale_x_reverse(limits = c(94.6, 94.3)) +
  #scale_x_continuous()+
  scale_y_continuous() +
  #scale_x_continuous(limits = c(94.6, 94.55)) +
  #coord_cartesian(xlim=c(-20, 80))+
  coord_flip()

ggplotly(pools_alttime)

CO2 calculations by date

alldata_w_time %>% filter(compound == "phytane") %>% select(Date.Ma, eps_p_lipid, d13C_pool_avg, true_d13C_pred, CO2_aq_conc, pCO2) %>% kable() 

Date.Ma	eps_p_lipid	d13C_pool_avg	true_d13C_pred	CO2_aq_conc	pCO2
94.45489	17.93042	-6.505161	-30.40516	24.04670	994.4872
94.30791	17.92757	-6.349668	-30.24967	24.03699	994.0855
94.34293	17.91180	-5.490317	-29.39032	23.98351	991.8741
94.22115	17.88922	-4.257208	-28.15721	23.90736	988.7246
94.37841	17.87707	-3.592498	-27.49250	23.86659	987.0385
94.17981	17.90731	-5.245692	-29.14569	23.96835	991.2470
94.04433	17.87767	-3.625275	-27.52528	23.86860	987.1214
94.56738	17.91410	-5.615702	-29.51570	23.99130	992.1959
94.52951	17.91765	-5.809358	-29.70936	24.00333	992.6935
93.63407	17.92682	-6.309279	-30.20928	24.03447	993.9813
93.91610	17.90231	-4.972727	-28.87273	23.95147	990.5487
93.66267	17.92183	-6.037430	-29.93743	24.01752	993.2804
94.27372	17.89758	-4.714443	-28.61444	23.93552	989.8892
93.85049	17.89374	-4.504509	-28.40451	23.92258	989.3540
94.10260	17.88144	-3.831769	-27.73177	23.88124	987.6445
94.01944	17.87633	-3.551591	-27.45159	23.86409	986.9350
93.78417	17.90732	-5.246106	-29.14611	23.96838	991.2481
93.73495	17.92271	-6.085177	-29.98518	24.02049	993.4034
93.68296	17.91013	-5.399526	-29.29953	23.97788	991.6412
93.75582	17.91882	-5.872967	-29.77297	24.00728	992.8571
93.61487	17.91623	-5.731822	-29.63182	23.99851	992.4942
93.57814	17.90877	-5.324885	-29.22488	23.97326	991.4499
93.81802	17.91766	-5.810122	-29.71012	24.00338	992.6954
94.13923	17.89399	-4.517941	-28.41794	23.92341	989.3882

eps_p

#using my eps_p (lipd rather than caco3)
eps_p <- final_data %>% 
  filter(compound %in% c( "phytane" , "C19", "C21"), true_d13C_pred_in_range == TRUE)%>% 
  iso_plot_data(
    x = depth, 
    y = eps_p_lipid 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "eps_p") +
  #scale_x_reverse() +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(eps_p)

pCO2

#using my eps_p (lipd rather than caco3)
pCO2 <- alldata_w_time %>% 
  filter(compound %in% c( "phytane" , "C19", "C21"), true_d13C_pred_in_range == TRUE)%>% 
  iso_plot_data(
    x = Date.Ma, 
    y = pCO2 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "pCO2") +
  scale_x_reverse() +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 
ggplotly(pCO2)

pCO22 <- alldata_w_time %>% 
  filter(compound %in% c( "C21"), true_d13C_pred_in_range == TRUE)%>% 
  ggplot()+
  aes(
    x = Date.Ma, 
    y = pCO2 
    #y_error = true_d13C_pred_se ,
    #shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_smooth() +
  geom_point() +
  #theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "pCO2") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(-100, 120))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 


alkane <- alldata_w_time %>% 
  filter(compound %in% c( "C31", "C33"), true_d13C_pred_in_range == TRUE, prep == "adduct")%>%
  ggplot()+
  aes(
    x = Date.Ma, 
    y = d13C_pool_avg, 
    #y_error = true_d13C_pred_se ,
    shape = compound 
    #color = prep
  ) + 
  # additional features beyond the default plot
  facet_grid(~compound, scales = "free") + 
  geom_smooth() +
  geom_point() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5)) +
  labs(y = "d13C") +
  scale_x_reverse(limits = c(94.5, 94.34)) +
  #scale_x_continuous(limits = c(94.6, 94.3))+
  #scale_y_continuous(limits = c(-33, -26)) +
  coord_flip() 

ggplotly(pCO2)

ggplotly(alkane)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'

grid.arrange( alkane, pCO22, ncol = 2)
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#> Warning: Removed 20 rows containing non-finite values (stat_smooth).
#> Warning: Removed 20 rows containing missing values (geom_point).